      subroutine parset
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx
      use blcom_com_M
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer , dimension(8) :: ijkv
      integer , dimension(27) :: ijkc, ijkstep
      integer :: n, np, kr, iref, istp, jstep, kstep, icount, l, kx, ky, kz, &
         ijkreg, iter, itmax_map, nbin, ij, ix, iz
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5, shaperx,shapery, bbb, vvv
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double), dimension(8) :: regvol
      real(double), dimension(27) :: wghtc
      real(double) :: nuxp, nuyp, nuzp, riwid, rjwid, rkwid, volmin, wsnc, &
         rnpcx, rnpcy, rnpcz, ws, xi, eta, zta, epsreg, wsxi, wseta, wsnu, &
         pvolume, flvx, flvy, flvz, zone, ztwo, variation, segno, harris, fx, &
         fi, vmag, th, cosfi, sinfi, shaper, bxp, byp, bzp, rbmag, udotb, &
         qptotl, qion, tionx, tiony, tionz, qele, telex, teley, telez, pion, &
         pele, wsq, wsqc, xc, zc,bfldx,bfldy,bfldz, t_shape,efldx,efldy,efldz
      real(double) :: gamma, px1, py1, pz1, gamma1, flvxp,flvyp,flvzp
      logical :: expnd, nomore
      character :: name*80
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      real(double) , external :: ranf
!-----------------------------------------------
!
!
!

      iplost = 0
      riwid = 1./float(iwid)
      rjwid = 1./float(jwid)
      rkwid = 1./float(kwid)
!
      call celstep (iwid, jwid, kwid, ijkstep)
!
!
!     set up linked list
!
      iphead(ijkcell(:ncells)) = 0
!
      iphead(1) = 1
      do np = 1, npart
         pxi(np) = 1.
         peta(np) = 1.
         pzta(np) = 1.
!
         link(np) = np + 1
      end do
      link(npart) = 0
!
!
!     set all particle variables for np=0
!
      link(0) = 0
      px(0) = 0.0
      py(0) = 0.0
      pz(0) = 0.0
      pmu(0) = 0.0
      pxi(0) = 0.0
      peta(0) = 0.0
      pzta(0) = 0.0
      up(0) = 0.0
      vp(0) = 0.0
      wp(0) = 0.0
      ico(0) = 0
      qpar(0) = 0.0
!
!     uniform distribution of particles throughout the volume
!     with random velocities
!
      npsampl = 0
!
!dir$ ivdep
      do n = 1, nvtx
         ijk = ijkvtx(n)
         a11(ijk) = jx0(ijk)/vvol(ijk)
         a12(ijk) = jy0(ijk)/vvol(ijk)
         a13(ijk) = jz0(ijk)/vvol(ijk)
      end do
!
if (initial_equilibrium.eq.6) then
open(15,file=nome_file_input,form='formatted',status='unknown')

do ij = 1, ibar*kbar
!
   read(15,*)ix,iz,bbb,vvv
!   write(*,*)'read  i=',ix,'   j=',iz,'   By=',bbb
   i=ix+1
   k=iz+1
!
   do j = 2, jbp1
!
      ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
        a21(ijk) = vvv
!
!
   end do
end do
close(15)
endif
!
      do kr = 1, nrg
!
!     check convexity of region
!
         iref = 8*(kr - 1)
         istp = 1
         jstep = 2
         kstep = 4
!
         call triple (xvi, yvi, zvi, iref + 1, 1, 3, 4, regvol(1))
         call triple (xvi, yvi, zvi, iref + 2, 1, -1, 4, regvol(2))
         call triple (xvi, yvi, zvi, iref + 3, 1, -1, 4, regvol(3))
         call triple (xvi, yvi, zvi, iref + 4, -3, -1, 4, regvol(4))
         call triple (xvi, yvi, zvi, iref + 5, 1, 3, -4, regvol(5))
         call triple (xvi, yvi, zvi, iref + 6, 1, -1, -4, regvol(6))
         call triple (xvi, yvi, zvi, iref + 7, 1, -1, -4, regvol(7))
         call triple (xvi, yvi, zvi, iref + 8, -3, -1, -4, regvol(8))
!
         regvol(5:8) = -regvol(5:8)
!
         volmin = 0.0
         do l = 1, 8
            volmin = amin1(regvol(l),volmin)
         end do
!
         if (volmin < 0.0) then
            write (*, *) 'region', nrg, ' is not convex'
            cycle
         endif
!
         wsnc = 1./(npcelx(kr)*npcely(kr)*npcelz(kr))
         rnpcx = 1./float(npcelx(kr))
         rnpcy = 1./float(npcely(kr))
         rnpcz = 1./float(npcelz(kr))
!
!     permute list of region vertices to lexicographic order
!
         ws = xvi(3,kr)
         xvi(3,kr) = xvi(4,kr)
         xvi(4,kr) = ws
!
         ws = yvi(3,kr)
         yvi(3,kr) = yvi(4,kr)
         yvi(4,kr) = ws
!
         ws = zvi(3,kr)
         zvi(3,kr) = zvi(4,kr)
         zvi(4,kr) = ws
!
         ws = xvi(7,kr)
         xvi(7,kr) = xvi(8,kr)
         xvi(8,kr) = ws
!
         ws = yvi(7,kr)
         yvi(7,kr) = yvi(8,kr)
         yvi(8,kr) = ws
!
         ws = zvi(7,kr)
         zvi(7,kr) = zvi(8,kr)
         zvi(8,kr) = ws
!
!
         do kx = 1, npcelx(kr)
            do ky = 1, npcely(kr)
               do kz = 1, npcelz(kr)
!
                  xi = (0.5 + (kx - 1))*rnpcx
                  eta = (0.5 + (ky - 1))*rnpcy
                  zta = (0.5 + (kz - 1))*rnpcz
!
                  call watec (xi, eta, zta, wghtc)
!
!
                  call weights (xi, eta, zta, wght)
!
                  ijkreg = 8*(kr - 1) + 1
                  istp = 1
                  jstep = 2
                  kstep = 4
                  epsreg = 1.E-4
!
                  do n = 1, ncells
!
                     ijk = ijkcell(n)
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!
                     xptilde(ijk) = wght(1)*x(ipjk) + (wght(2)*x(ipjpk)+(wght(3&
                        )*x(ijpk)+(wght(4)*x(ijk)+(wght(5)*x(ipjkp)+(wght(6)*x(&
                        ipjpkp)+(wght(7)*x(ijpkp)+wght(8)*x(ijkp)))))))
!
                     yptilde(ijk) = wght(1)*y(ipjk) + (wght(2)*y(ipjpk)+(wght(3&
                        )*y(ijpk)+(wght(4)*y(ijk)+(wght(5)*y(ipjkp)+(wght(6)*y(&
                        ipjpkp)+(wght(7)*y(ijpkp)+wght(8)*y(ijkp)))))))
!
                     zptilde(ijk) = wght(1)*z(ipjk) + (wght(2)*z(ipjpk)+(wght(3&
                        )*z(ijpk)+(wght(4)*z(ijk)+(wght(5)*z(ipjkp)+(wght(6)*z(&
                        ipjpkp)+(wght(7)*z(ijpkp)+wght(8)*z(ijkp)))))))
!
                  end do
!
!
!     check whether particle is in region or not
!
                  do n = 1, ncells
                     np = iphead(1)
                     if (np == 0) then
                     write(*,*)'Hello Dave, you did not put enough particles. I am stopping', &
                                ' I suggest you take a space walk and change the number of particles in corgan'
                     stop
                     end if
!
                     ijk = ijkcell(n)
                     call map3d (ijkreg, istp, jstep, kstep, epsreg, iter, &
                        itmax_map, xvi, yvi, zvi, xptilde(ijk), yptilde(ijk), &
                        zptilde(ijk), wsxi, wseta, wsnu)
                     if (wsxi*(1. - wsxi)>=0.0 .and. wseta*(1.-wseta)>=0.0&
                         .and. wsnu*(1.-wsnu)>=0.0) then
!
                        iphead(1) = link(np)
                        link(np) = iphead(ijk)
                        iphead(ijk) = np
!
                        k = 1 + int((ijk - 1)*rkwid)
                        j = 1 + int((ijk - 1 - (k - 1)*kwid)*rjwid)
                        i = 1 + int((ijk - 1 - (j - 1)*jwid - (k - 1)*kwid)*&
                           riwid)
!
!
                        pxi(np) = i + xi
                        peta(np) = j + eta
                        pzta(np) = k + zta
!
                        px(np) = xptilde(ijk)
                        py(np) = yptilde(ijk)
                        pz(np) = zptilde(ijk)
!
                     else
!
                        !write (*, *) &
                        !   'The input file has a problem in desinging the particle domain',wsxi,wseta,wsnu,np
                        !stop
!
                     endif
!
                  end do
!
!     calculate correct velocity distribution for 3 velocity components
!
                  call gauss3v (nsampl, vrms, proble, wate(1,1), wate(1,2), &
                     wate(1,3), wate(1,4))
!
!
!dir$ ivdep
                  do n = 1, ncells
!
                     np = iphead(ijkcell(n))
!
                     ijk = ijkcell(n)
!
!     use Sulsky's routine to calculate partial volumes
!
                     call parvol (kx, ky, kz, rnpcx, rnpcy, rnpcz, x, y, z, ijk&
                        , iwid, jwid, kwid, pvolume)
!
!
                     ipjk = ijk + iwid
                     ipjpk = ijk + iwid + jwid
                     ijpk = ijk + jwid
!
                     ijkp = ijk + kwid
                     ipjkp = ijk + iwid + kwid
                     ijpkp = ijk + jwid + kwid
                     ipjpkp = ijk + iwid + jwid + kwid
!
                     segno = sign(1.0,qom(icoi(kr)))
!
!       Particle fluid velocity interpolated from the
!       current
!

                    if(fluid_load) then

                     call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)

!                     flvx = wght(1)*a11(ipjk) + (wght(2)*a11(ipjpk)+(wght(3)*&
!                        a11(ijpk)+(wght(4)*a11(ijk)+(wght(5)*a11(ipjkp)+(wght(6&
!                        )*a11(ipjpkp)+(wght(7)*a11(ijpkp)+wght(8)*a11(ijkp)))))&
!                        ))
!
!                     flvy = wght(1)*a12(ipjk) + (wght(2)*a12(ipjpk)+(wght(3)*&
!                        a12(ijpk)+(wght(4)*a12(ijk)+(wght(5)*a12(ipjkp)+(wght(6&
!                        )*a12(ipjpkp)+(wght(7)*a12(ijpkp)+wght(8)*a12(ijkp)))))&
!                        ))
!
!                     flvz = wght(1)*a13(ipjk) + (wght(2)*a13(ipjpk)+(wght(3)*&
!                        a13(ijpk)+(wght(4)*a13(ijk)+(wght(5)*a13(ipjkp)+(wght(6&
!                        )*a13(ipjpkp)+(wght(7)*a13(ijpkp)+wght(8)*a13(ijkp)))))&
!                        ))
!
!                    ASSUMES all current carried by electrons
!
!                     if(segno.lt.0.0) then
!                     flvx = -flvx/rhr(kr)/harris/2
!                     flvy = -flvy/rhr(kr)/harris/2
!                     flvz = -flvz/rhr(kr)/harris/2
!                     else
!                     flvx = flvx/rhr(kr)/harris/2
!                     flvy = flvy/rhr(kr)/harris/2
!                     flvz = flvz/rhr(kr)/harris/2
!                     endif
!
                     else
!
                     call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
!
                     end if
!
                     ico(np) = icoi(kr)
!
                     qpar(np) = pvolume*rhr(kr)*segno*harris
!
                     fi = acos(2.*ranf() - 1.)
                     nbin = int(ranf()*float(nsampl - 1) + 1.)
                     vmag = vrms(nbin)*sqrt(2.)*t_shape
                     th = 2.*pi*ranf()
                     cosfi = 2.*ranf() - 1.
                     sinfi = sqrt(1. - cosfi**2)
                     if (ranf() < 0.5) sinfi = -sinfi
!
 					if(relativity) then

						call add_relativ(clite,uvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np)), &
									     vvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np)), &
									     wvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np)), &
                                         flvx,flvy,flvz,flvxp,flvyp,flvzp)

 						px1 = siepx(kr)*vmag*cos(th)*sinfi

                        py1 = siepy(kr)*vmag*sin(th)*sinfi

                        pz1 = siepz(kr)*vmag*cosfi

                        gamma1 = sqrt(1.0+px1**2+py1**2+pz1**2)

                        call add_relativ(clite,px1/gamma1,py1/gamma1,pz1/gamma1, &
                                         flvxp,flvyp,flvzp,up(np),vp(np),wp(np))

!
!                        up(np) = (flvx + siepx(kr)*vmag*cos(th)*sinfi)/&
!                                 (1.0+flvx * siepx(kr)*vmag*cos(th)*sinfi)
!
!                        vp(np) = (flvy + siepy(kr)*vmag*sin(th)*sinfi)/&
!                                 (1.0+ flvy * siepy(kr)*vmag*sin(th)*sinfi)
!
!                        wp(np) = (flvz + siepz(kr)*vmag*cosfi)/&
!                                 (1.0+ flvz * siepz(kr)*vmag*cosfi)

					else


                        flvx=flvx+ uvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))
                        flvy=flvy+ vvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))
                        flvz=flvz+ wvip(kr)*cos(2.*pi/(kbar*dz)*mpert*pz(np))

						up(np) = flvx + siepx(kr)*vmag*cos(th)*sinfi

                        vp(np) = flvy + siepy(kr)*vmag*sin(th)*sinfi

                        wp(np) = flvz + siepz(kr)*vmag*cosfi



					end if

if (initial_equilibrium.eq.6) then
   vp(np)=vp(np)+a21(ijk)
end if

!	if(relativity) then
!           gamma = sqrt(1.0+up(np)**2+vp(np)**2+wp(np)**2)
!           up(np) = up(np)/gamma
!           vp(np) = vp(np)/gamma
!           wp(np) = wp(np)/gamma
!	end if
!
!     calculate particle magnetic moment
!
                     call watec (pxi(np), peta(np), pzta(np), wghtc)
                     i = int(pxi(np))
                     j = int(peta(np))
                     k = int(pzta(np))
!
                     call celdex (i, j, k, iwid, jwid, kwid, ijkc)
!
                     call triquad (ijkc, bxn, byn, bzn, wghtc, bxp, byp, bzp)
!
                     rbmag = 1./(sqrt(bxp**2 + byp**2 + bzp**2) + 1.E-20)
!
                     udotb = (up(np)*bxp+vp(np)*byp+wp(np)*bzp)*rbmag
!
                     pmu(np) = 0.5*(up(np)**2+vp(np)**2+wp(np)**2-udotb**2)*&
                        rbmag
!
!
                  end do
               end do
            end do
         end do
      end do
!
!     compute total particle charge and temperature by species
!
      qptotl = 0.0
      qion = 0.
      tionx = 0.
      tiony = 0.
      tionz = 0.
      qele = 0.
      telex = 0.
      teley = 0.
      telez = 0.
      pion = 0.
      pele = 0.
      do n = 1, ncells
         ijk = ijkcell(n)
         np = iphead(ijk)
         wsq = 0.0
 2001    continue
         if (np > 0) then
            wsq = wsq + qpar(np)
            qptotl = qptotl + qpar(np)
            segno = sign(1.0,qpar(np))
            if (ico(np) == 1) then
               qion = qion + qpar(np)
               tionx = tionx + qpar(np)*(up(np)-uvi(1)*segno)**2
               tiony = tiony + qpar(np)*(vp(np)-vvi(1)*segno)**2
               tionz = tionz + qpar(np)*(wp(np)-wvi(1)*segno)**2
               pion = pion + qpar(np)*vp(np)
            else
               qele = qele + qpar(np)
               telex = telex + qpar(np)*(up(np)+uvi(2)*segno)**2
               teley = teley + qpar(np)*(vp(np)+vvi(2)*segno)**2
               telez = telez + qpar(np)*(wp(np)+wvi(2)*segno)**2
               pele = pele + qpar(np)*vp(np)
            endif
            np = link(np)
            go to 2001
         endif
         wsqc = wsq/vol(ijk)
      end do
      telex = telex/qele/abs(qom(2))
      teley = teley/qele/abs(qom(2))
      telez = telez/qele/abs(qom(2))
      tionx = tionx/qion/abs(qom(1))
      tiony = tiony/qion/abs(qom(1))
      tionz = tionz/qion/abs(qom(1))
      pion = pion/qion
      pele = pele/qele
      write (*, *) 'parset: qptotl=', qptotl
      write (*, *) 'parset: elet x=', telex
      write (*, *) 'parset: elet y=', teley
      write (*, *) 'parset: elet z=', telez
      write (*, *) 'parset: iont x=', tionx
      write (*, *) 'parset: iont y=', tiony
      write (*, *) 'parset: iont z=', tionz
      write (*, *) 'parset: p ele y=', pele
      write (*, *) 'parset: p ion y=', pion
!
!
      return
!
      end subroutine parset

      subroutine add_relativ(clite,ux,uy,uz,vx,vy,vz,sx,sy,sz)
      USE vast_kind_param, ONLY:  double
      implicit none
      real(double) :: clite,ux,uy,uz,vx,vy,vz,sx,sy,sz
      real(double) :: uparx,upary,uparz,uperpx,uperpy,uperpz,vmod2,utmp,alfav,denom

      vmod2=(vx**2+vy**2+vz**2)

      utmp=(vx*ux+vy*uy+vz*uz)/(vmod2+1d-10)

      uparx=utmp*vx
      upary=utmp*vy
      uparz=utmp*vz

      uperpx=ux-uparx
      uperpy=uy-upary
      uperpz=uz-uparz

      alfav=sqrt(1.0-(vx**2+vy**2+vz**2)/clite**2)
      denom=1+(vx*ux+vy*uy+vz*uz)/clite**2

      sx=(vx+uparx+alfav*uperpx)/denom
      sy=(vy+upary+alfav*uperpy)/denom
      sz=(vz+uparz+alfav*uperpz)/denom

      end subroutine add_relativ
